Quad and Hami SR results
Mick Girdwood
06 December, 2023
library(metafor) # for meta-analysis
library(clubSandwich) # for covariance matrix and robust ci estimates
library(rms) # for fitting splines
library(ellipse)
library(ggpubr)
library(mgfunctions)
library(kableExtra)1. Quadriceps Strength - Within Person
# Generate the within person quad strength data
quad <- within_data %>%
rename(acl_graft_group = graft_group) %>%
filter(str_detect(measure, "quad")) %>% # take all quad based outcomes
filter(timepoint_mean >2, # remove pre-operative data
timepoint_mean < 120, # remove >10 year data
str_detect(graft, "contralateral|Contralateral", negate = TRUE)) %>% # remove any contralateral grafts
mutate(es_id = row_number()) %>% # effect size id number for use in random effects
mutate(measure_2 = case_when( # classify different outcomes into subgroups:
str_detect(measure, "mvic|hhd") ~ "Isometric",
str_detect(measure, "isk con 180|isk con 230|isk con 240|isk con 300") ~ "Fast Isokinetic",
str_detect(measure, "isk con 30|isk con 60|isk con 90|isk con 120") ~ "Slow Isokinetic"
)) %>%
filter(!is.na(measure_2))
# Reducing timepoint down to categories (3, 6, 12, 24, 48, 96 months post)
# not using this at the moment but keeping code here for now
quad_cat <- quad %>%
mutate(timepoint_cut = cut(timepoint_mean,
breaks = c(0, 4.5, 9, 18, 36, 72, Inf),
labels = c(3, 6, 12, 24, 48, 96)))
#group_by(cohort, timepoint_cut, group, measure_2) %>%
# if multiple timepoints allocated to same category, take the closest to the categorical timepoint
#slice(which.min(abs(timepoint_mean - as.numeric(as.character(timepoint_cut))))) %>%
#ungroup()
# separate subgroups of data based on contraction type
# Some studies report multiple outcomes in the same subgroup
# Case by case removal of data
# Generally:
# for fast if they have isk con 180 use that if possible, otherwise use speed closest to this.
# for slow use isk con 60, or 90
fastdata <- quad_cat %>% filter(measure_2 == "Fast Isokinetic") %>%
filter(!(cohort == "Bailey 2019" & measure == "quad isk con 300"),
!(cohort == "Drocco 2017" & measure == "quad isk con 240"),
!(cohort == "Kyritsis 2016" & measure == "quad isk con 300"),
!(cohort == "Laudner 2015" & measure == "quad isk con 300"),
!(cohort == "Tourville 2014" & measure == "quad isk con 300"),
!(cohort == "Welling 2020" & measure == "quad isk con 300"),
!(cohort == "Welling 2018b" & measure == "quad isk con 300")
)
slowdata <- quad_cat %>% filter(measure_2 == "Slow Isokinetic") %>%
filter(!(cohort == "Li 2022" & measure == "quad isk con 30"),
!(cohort == "Siney 2010" & measure == "quad isk con 90"),
!(cohort == "Zult 2017" & measure == "quad isk con 120"),
!(cohort == "Pamukoff 2018" & measure == "quad isk con 120"))
isodata <- quad_cat %>% filter(measure_2 == "Isometric") %>%
filter(!(cohort == "Karanikas 2005" & measure == "quad mvic 0"),
!(cohort == "Labanca 2018" & measure == "quad mvic 30"),
!(cohort == "Labanca 2016" & measure == "quad mvic 30"),
!(cohort == "Wongcharoenwatana 2019" & measure == "quad mvic 90 hhd prone"))Specify covariance matrix
Need to create a full covariance matrix. Using clubSandwhich package to help Assuming a correlation of 0.85 between timepoints.
fastV <- impute_covariance_matrix(vi = fastdata$vi,
cluster = fastdata$cohort, # cluster is cohort (not study)
ti = fastdata$timepoint_mean, # timepoint
ar1 = 0.85, # auto-correlation between timepoints
check_PD = TRUE, # check positive definite afterwards
smooth_vi = TRUE,
return_list = FALSE) # return the full matrix
slowV <- impute_covariance_matrix(vi = slowdata$vi,
cluster = slowdata$cohort,
ti = slowdata$timepoint_mean,
ar1 = 0.85,
smooth_vi = TRUE,
return_list = FALSE)
isoV <- impute_covariance_matrix(vi = isodata$vi,
cluster = isodata$cohort,
ti = isodata$timepoint_mean,
ar1 = 0.85,
check_PD = TRUE,
smooth_vi = TRUE,
return_list = FALSE)Model Selection
Data are now ready for meta-analysis
Need to now decide on how to fit models. Several different structures were trialed in piloting.
Based on piloting best results were achieved with fitting random effects for timepoints, nested within cohorts. Trialled fitting extra random effects e.g. effect size id (es_id), as well as separate models where study group was also a separate random effect, however resulting models were too complex, with indistinguishable random effects and overall a simpler model chosen.
Decision making here revolves around how to fit the timepoint predictor - i.e. what sort of relationship is present between yi and timepoint Different models are generated an then using fit statistics (AIC, BIC, AIcc), visual inspection of fit and expected fit based on knowledge to decide on best fit.
5 different shapes of fit tried: - linear - log - polynomial - 3 knot restricted cubic spline - 4 knot restricted cubic spline
# Create 'empty' models without timepoint used as a moderator variable.
# Will use these models to then build different versions on
# Using a continuous time auto-regressive variance structure (CAR)
slow_mv_empty <- rma.mv(yi, slowV,
data = slowdata,
random = list(~ timepoint_mean|cohort),
struct = c("CAR"))
fast_mv_empty <- rma.mv(yi, fastV,
data = fastdata,
random = list(~ timepoint_mean|cohort),
struct = c("CAR"))
iso_mv_empty <- rma.mv(yi, isoV,
data = isodata,
random = list(~ timepoint_mean|cohort),
struct = c("CAR"))Slow Isokinetic Quadriceps
4 knot spline best fit
mod_selection(slow_mv_empty)## [[1]]
## mod logLik. deviance. AIC. BIC. AICc.
## 1 Linear 150.1557 -300.3114 -292.3114 -278.2901 -292.1454
## 2 Log 190.8914 -381.7827 -373.7827 -359.7614 -373.6168
## 3 Poly (2) 165.4130 -330.8259 -320.8259 -303.3196 -320.5749
## 4 3 knot RCS 190.4994 -380.9989 -370.9989 -353.4926 -370.7478
## 5 4 knot RCS 225.2438 -450.4876 -438.4876 -417.5046 -438.1332
##
## [[2]]
Fast Isokinetic Quadriceps
4 knot spline best fit
mod_selection(fast_mv_empty)## [[1]]
## mod logLik. deviance. AIC. BIC. AICc.
## 1 Linear 83.57804 -167.1561 -159.1561 -148.8162 -158.7260
## 2 Log 92.82967 -185.6593 -177.6593 -167.3195 -177.2292
## 3 Poly (2) 88.06118 -176.1224 -166.1224 -153.2488 -165.4630
## 4 3 knot RCS 92.70016 -185.4003 -175.4003 -162.5268 -174.7410
## 5 4 knot RCS 98.19512 -196.3902 -184.3902 -169.0041 -183.4464
##
## [[2]]
Isometric Quadriceps
Log is best fitting, however 4 knot spline appears more realistic, with marginally worse fit.
mod_selection(iso_mv_empty)## [[1]]
## mod logLik. deviance. AIC. BIC. AICc.
## 1 Linear 26.12579 -52.25159 -44.25159 -35.94144 -43.51085
## 2 Log 32.34395 -64.68791 -56.68791 -48.37776 -55.94717
## 3 Poly (2) 28.34331 -56.68662 -46.68662 -36.38441 -45.53277
## 4 3 knot RCS 30.17380 -60.34761 -50.34761 -40.04539 -49.19376
## 5 4 knot RCS 32.94939 -65.89877 -53.89877 -41.64047 -52.21877
##
## [[2]]
Final models
Slow Isokinetic Quadriceps
slow_mv <- rma.mv(yi, slowV,
mods = ~rcs(timepoint_mean, 4),
data = slowdata,
random = list(~ timepoint_mean|cohort),
struct = c("CAR"))
robust(slow_mv, cluster = cohort, clubSandwich = TRUE)##
## Multivariate Meta-Analysis Model (k = 248; method: REML)
##
## Variance Components:
##
## outer factor: cohort (nlvls = 165)
## inner factor: timepoint_mean (nlvls = 106)
##
## estim sqrt fixed
## tau^2 0.0108 0.1038 no
## rho 0.9209 no
##
## Test for Residual Heterogeneity:
## QE(df = 244) = 8519.4758, p-val < .0001
##
## Number of estimates: 248
## Number of clusters: 165
## Estimates per cluster: 1-7 (mean: 1.50, median: 1)
##
## Test of Moderators (coefficients 2:4):¹
## F(df1 = 3, df2 = 31.35) = 82.1206, p-val < .0001
##
## Model Results:
##
## estimate se¹ tval¹ df¹
## intrcpt -0.7034 0.0399 -17.6451 36.59
## rcs(timepoint_mean, 4)timepoint_mean 0.0749 0.0063 11.9131 39.34
## rcs(timepoint_mean, 4)timepoint_mean' -2.2454 0.2212 -10.1528 50.52
## rcs(timepoint_mean, 4)timepoint_mean'' 3.3612 0.3348 10.0408 51.58
## pval¹ ci.lb¹ ci.ub¹
## intrcpt <.0001 -0.7842 -0.6226 ***
## rcs(timepoint_mean, 4)timepoint_mean <.0001 0.0622 0.0876 ***
## rcs(timepoint_mean, 4)timepoint_mean' <.0001 -2.6895 -1.8013 ***
## rcs(timepoint_mean, 4)timepoint_mean'' <.0001 2.6893 4.0330 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## 1) results based on cluster-robust inference (var-cov estimator: CR2,
## approx t/F-tests and confidence intervals, df: Satterthwaite approx)profile(slow_mv)## Profiling tau2 = 1
## Profiling rho = 1
Plots - normal and log scale
mv_plotdetails(slow_mv, include_pi = TRUE)
mv_plotdetails(slow_mv, logscale = TRUE, showgraft = FALSE)
predict_details(slow_mv) %>% kbl() %>% kable_styling(position = "left", full_width = FALSE)| 1 year | 2 years | 5 years | Zero crossing | Last Data Point |
|---|---|---|---|---|
| 84.9 (83.5 to 86.3) | 87.4 (85.7 to 89.1) | 91.3 (89.2 to 93.4) | 8.8 years | 9 years |
Fast Isokinetic Quadriceps
fast_mv <- rma.mv(yi, fastV,
mods = ~rcs(timepoint_mean, 4),
data = fastdata,
random = list(~ timepoint_mean|cohort),
struct = c("CAR"))
robust(fast_mv, cluster = cohort, clubSandwich = TRUE)##
## Multivariate Meta-Analysis Model (k = 100; method: REML)
##
## Variance Components:
##
## outer factor: cohort (nlvls = 66)
## inner factor: timepoint_mean (nlvls = 43)
##
## estim sqrt fixed
## tau^2 0.0091 0.0952 no
## rho 0.9298 no
##
## Test for Residual Heterogeneity:
## QE(df = 96) = 2868.2638, p-val < .0001
##
## Number of estimates: 100
## Number of clusters: 66
## Estimates per cluster: 1-6 (mean: 1.52, median: 1)
##
## Test of Moderators (coefficients 2:4):¹
## F(df1 = 3, df2 = 16.1) = 15.5012, p-val < .0001
##
## Model Results:
##
## estimate se¹ tval¹ df¹
## intrcpt -0.5051 0.0595 -8.4955 14.86
## rcs(timepoint_mean, 4)timepoint_mean 0.0489 0.0097 5.0487 15.72
## rcs(timepoint_mean, 4)timepoint_mean' -1.1303 0.2765 -4.0879 20.48
## rcs(timepoint_mean, 4)timepoint_mean'' 1.6862 0.4200 4.0151 20.99
## pval¹ ci.lb¹ ci.ub¹
## intrcpt <.0001 -0.6319 -0.3783 ***
## rcs(timepoint_mean, 4)timepoint_mean 0.0001 0.0283 0.0695 ***
## rcs(timepoint_mean, 4)timepoint_mean' 0.0005 -1.7062 -0.5544 ***
## rcs(timepoint_mean, 4)timepoint_mean'' 0.0006 0.8128 2.5595 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## 1) results based on cluster-robust inference (var-cov estimator: CR2,
## approx t/F-tests and confidence intervals, df: Satterthwaite approx)profile(fast_mv)## Profiling tau2 = 1
## Profiling rho = 1
Plots - normal and log scale
mv_plotdetails(fast_mv, include_pi = TRUE)
mv_plotdetails(fast_mv, logscale = TRUE, showgraft = FALSE)
predict_details(fast_mv) %>% kbl() %>% kable_styling(position = "left", full_width = FALSE)| 1 year | 2 years | 5 years | Zero crossing | Last Data Point |
|---|---|---|---|---|
| 86.4 (84.4 to 88.5) | 88.7 (86.5 to 91) | 89.9 (84 to 96.2) | 6.5 years | 5.1 years |
Isometric Quadriceps
iso_mv <- rma.mv(yi, isoV,
mods = ~rcs(timepoint_mean, 4),
data = isodata,
random = list(~ timepoint_mean|cohort),
struct = c("CAR"))
robust(iso_mv, cluster = cohort, clubSandwich = TRUE)##
## Multivariate Meta-Analysis Model (k = 61; method: REML)
##
## Variance Components:
##
## outer factor: cohort (nlvls = 45)
## inner factor: timepoint_mean (nlvls = 58)
##
## estim sqrt fixed
## tau^2 0.0179 0.1337 no
## rho 0.0000 no
##
## Test for Residual Heterogeneity:
## QE(df = 57) = 3768.7603, p-val < .0001
##
## Number of estimates: 61
## Number of clusters: 45
## Estimates per cluster: 1-3 (mean: 1.36, median: 1)
##
## Test of Moderators (coefficients 2:4):¹
## F(df1 = 3, df2 = 8.97) = 14.7038, p-val = 0.0008
##
## Model Results:
##
## estimate se¹ tval¹ df¹
## intrcpt -0.5719 0.0825 -6.9351 13.57
## rcs(timepoint_mean, 4)timepoint_mean 0.0560 0.0127 4.4242 18.53
## rcs(timepoint_mean, 4)timepoint_mean' -0.8114 0.2265 -3.5820 24.57
## rcs(timepoint_mean, 4)timepoint_mean'' 1.0613 0.3010 3.5261 25.03
## pval¹ ci.lb¹ ci.ub¹
## intrcpt <.0001 -0.7493 -0.3945 ***
## rcs(timepoint_mean, 4)timepoint_mean 0.0003 0.0295 0.0826 ***
## rcs(timepoint_mean, 4)timepoint_mean' 0.0015 -1.2784 -0.3445 **
## rcs(timepoint_mean, 4)timepoint_mean'' 0.0017 0.4415 1.6812 **
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## 1) results based on cluster-robust inference (var-cov estimator: CR2,
## approx t/F-tests and confidence intervals, df: Satterthwaite approx)profile(iso_mv)## Profiling tau2 = 1
## Profiling rho = 1
Plots - normal and log scale
mv_plotdetails(iso_mv, include_pi = TRUE)
mv_plotdetails(iso_mv, logscale = TRUE, showgraft = FALSE)
predict_details(iso_mv) %>% kbl() %>% kable_styling(position = "left", full_width = FALSE)| 1 year | 2 years | 5 years | Zero crossing | Last Data Point |
|---|---|---|---|---|
| 90.9 (85.9 to 96.2) | 92.2 (87.5 to 97) | 99.2 (87 to 113.1) | 3.6 years | 5.4 years |
##########3
2. Hamstring Strength - Within Person
# Data for Hamstring
hs <- within_data %>%
rename(acl_graft_group = graft_group) %>%
filter(str_detect(measure, "hs")) %>%
filter(timepoint_mean >2, # remove pre-operative data
timepoint_mean < 120, # remove >10 year data
str_detect(graft, "contralateral|Contralateral", negate = TRUE)) %>% # remove any contralateral grafts
mutate(es_id = row_number()) %>%
mutate(measure_2 = case_when(
str_detect(measure, "mvic|hhd") ~ "Isometric",
str_detect(measure, "isk con 180|isk con 230|isk con 240|isk con 300") ~ "Fast Isokinetic",
str_detect(measure, "isk con 30|isk con 60|isk con 90|isk con 120") ~ "Slow Isokinetic"
)) %>%
filter(!is.na(measure_2))
# Reducing timepoint down to categories (3, 6, 12, 24, 48, 96 months post)
# not using this at the moment but keeping code here for now
hs_cat <- hs %>%
mutate(timepoint_cut = cut(timepoint_mean,
breaks = c(0, 4.5, 9, 18, 36, 72, Inf),
labels = c(3, 6, 12, 24, 48, 96)))
#group_by(cohort, timepoint_cut, measure_2) %>%
# if multiple timepoints allocated to same category, take the closest to the categorical timepoint
#slice(which.min(abs(timepoint_mean - as.numeric(as.character(timepoint_cut))))) %>%
#ungroup()
# separate subgroups of data based on contraction type
# Some studies report multiple outcomes in the same subgroup
# Case by case removal of data
# Generally:
# for fast if they have isk con 180 use that if possible, otherwise use speed closest to this.
# for slow use isk con 60, or 90
hs_fastdata <- hs_cat %>% filter(measure_2 == "Fast Isokinetic") %>%
filter(!(cohort == "Drocco 2017" & measure == "hs isk con 240"),
!(cohort == "Kyritsis 2016" & measure == "hs isk con 300"),
!(cohort == "Laudner 2015" & measure == "hs isk con 300"),
!(cohort == "Tourville 2014" & measure == "hs isk con 300"),
!(cohort == "Welling 2020" & measure == "hs isk con 300"),
!(cohort == "Welling 2018b" & measure == "hs isk con 300"))
hs_slowdata <- hs_cat %>% filter(measure_2 == "Slow Isokinetic") %>%
filter(!(cohort == "Lee 2019" & measure == "hs isk con 60 - deep"),
!(cohort == "Kim 2011" & measure == "hs isk con 90 hyperflex"),
!(cohort == "Siney 2010" & measure == "hs isk con 90"),
!(cohort == "McRae 2013" & measure == "hs isk con 60 supine"),
!(cohort == "Li 2022" & measure == "hs isk con 30"),
!(cohort == "Jiang 2012" & measure == "hs isk con 120"),
!(cohort == "Zult 2017" & measure == "hs isk con 120")
)
hs_isodata <- hs_cat %>% filter(measure_2 == "Isometric") %>%
filter(!(cohort == "Hu 2020" & measure == "hs mvic 30 hhd"),
!(cohort == "Hu 2020" & measure == "hs mvic 60 hhd"),
!(cohort == "Ardern 2010" & measure == "hs mvic 30"),
!(cohort == "Ardern 2010" & measure == "hs mvic 105"))Specify covariance matrix
Need to create a full covariance matrix. Using clubSandwhich package to help Assuming a correlation of 0.85 between timepoints.
hs_fastV <- impute_covariance_matrix(vi = hs_fastdata$vi,
cluster = hs_fastdata$cohort,
ti = hs_fastdata$timepoint_mean,
ar1 = 0.85,
check_PD = TRUE,
smooth_vi = TRUE,
return_list = FALSE)
hs_slowV <- impute_covariance_matrix(vi = hs_slowdata$vi,
cluster = hs_slowdata$cohort,
ti = hs_slowdata$timepoint_mean,
ar1 = 0.85,
smooth_vi = TRUE,
return_list = FALSE)
hs_isoV <- impute_covariance_matrix(vi = hs_isodata$vi,
cluster = hs_isodata$cohort,
ti = hs_isodata$timepoint_mean,
ar1 = 0.85,
check_PD = TRUE,
smooth_vi = TRUE,
return_list = FALSE)Model Selection
# Create 'empty' models without timepoint used as a moderator variable.
# Will use these models to then build different versions on
# Using a continuous time auto-regressive variance structure (CAR)
hs_fast_mv_empty <- rma.mv(yi, hs_fastV,
data = hs_fastdata,
random = list(~ timepoint_mean|cohort),
struct = c("CAR"))
hs_slow_mv_empty <- rma.mv(yi, hs_slowV,
data = hs_slowdata,
random = list(~ timepoint_mean|cohort),
struct = c("CAR"))
hs_iso_mv_empty <- rma.mv(yi, hs_isoV,
data = hs_isodata,
random = list(~ timepoint_mean|cohort),
struct = c("CAR"))Slow Isokinetic Hamstrings
4 knot spline best fit
mod_selection(hs_slow_mv_empty)## [[1]]
## mod logLik. deviance. AIC. BIC. AICc.
## 1 Linear 160.7531 -321.5061 -313.5061 -300.5817 -313.2864
## 2 Log 172.0767 -344.1534 -336.1534 -323.2290 -335.9337
## 3 Poly (2) 162.9579 -325.9158 -315.9158 -299.7871 -315.5825
## 4 3 knot RCS 170.9374 -341.8748 -331.8748 -315.7460 -331.5414
## 5 4 knot RCS 193.1509 -386.3018 -374.3018 -354.9796 -373.8299
##
## [[2]]
Fast Isokinetic Hamstrings
4 knot spline best fit
mod_selection(hs_fast_mv_empty)## [[1]]
## mod logLik. deviance. AIC. BIC. AICc.
## 1 Linear 70.06780 -140.1356 -132.1356 -122.7604 -131.5800
## 2 Log 72.93478 -145.8696 -137.8696 -128.4943 -137.3140
## 3 Poly (2) 70.64831 -141.2966 -131.2966 -119.6430 -130.4395
## 4 3 knot RCS 71.29610 -142.5922 -132.5922 -120.9385 -131.7351
## 5 4 knot RCS 75.16984 -150.3397 -138.3397 -124.4348 -137.1044
##
## [[2]]
Isometric hamstrings
Log is best fit 4 knot is overfit.
mod_selection(hs_iso_mv_empty)## [[1]]
## mod logLik. deviance. AIC. BIC. AICc.
## 1 Linear 11.84160 -23.68320 -15.68320 -9.349126 -14.39288
## 2 Log 12.53822 -25.07644 -17.07644 -10.742363 -15.78612
## 3 Poly (2) 11.13393 -22.26785 -12.26785 -4.491112 -10.19889
## 4 3 knot RCS 11.78507 -23.57014 -13.57014 -5.793402 -11.50118
## 5 4 knot RCS 16.69763 -33.39527 -21.39527 -12.237107 -18.28416
##
## [[2]]
Final models
Slow Isokinetic Hamstrings
hs_slow_mv <- rma.mv(yi, hs_slowV,
mods = ~rcs(timepoint_mean, 4),
data = hs_slowdata,
random = list(~ timepoint_mean|cohort),
struct = c("CAR"))
robust(hs_slow_mv, cluster = cohort, clubSandwich = TRUE)##
## Multivariate Meta-Analysis Model (k = 189; method: REML)
##
## Variance Components:
##
## outer factor: cohort (nlvls = 127)
## inner factor: timepoint_mean (nlvls = 90)
##
## estim sqrt fixed
## tau^2 0.0071 0.0844 no
## rho 0.8207 no
##
## Test for Residual Heterogeneity:
## QE(df = 185) = 5528.5858, p-val < .0001
##
## Number of estimates: 189
## Number of clusters: 127
## Estimates per cluster: 1-7 (mean: 1.49, median: 1)
##
## Test of Moderators (coefficients 2:4):¹
## F(df1 = 3, df2 = 24.88) = 11.5255, p-val < .0001
##
## Model Results:
##
## estimate se¹ tval¹ df¹
## intrcpt -0.4196 0.0569 -7.3787 27.59
## rcs(timepoint_mean, 4)timepoint_mean 0.0529 0.0089 5.9512 32.57
## rcs(timepoint_mean, 4)timepoint_mean' -1.7377 0.2982 -5.8268 40.88
## rcs(timepoint_mean, 4)timepoint_mean'' 2.6183 0.4501 5.8176 41.63
## pval¹ ci.lb¹ ci.ub¹
## intrcpt <.0001 -0.5361 -0.3030 ***
## rcs(timepoint_mean, 4)timepoint_mean <.0001 0.0348 0.0709 ***
## rcs(timepoint_mean, 4)timepoint_mean' <.0001 -2.3400 -1.1354 ***
## rcs(timepoint_mean, 4)timepoint_mean'' <.0001 1.7098 3.5269 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## 1) results based on cluster-robust inference (var-cov estimator: CR2,
## approx t/F-tests and confidence intervals, df: Satterthwaite approx)profile(hs_slow_mv)## Profiling tau2 = 1
## Profiling rho = 1
Plots - normal and log scale
mv_plotdetails(hs_slow_mv, include_pi = TRUE)
mv_plotdetails(hs_slow_mv, logscale = TRUE, showgraft = FALSE)
predict_details(hs_slow_mv) %>% kbl() %>% kable_styling(position = "left", full_width = FALSE)| 1 year | 2 years | 5 years | Zero crossing | Last Data Point |
|---|---|---|---|---|
| 93.8 (92.5 to 95.1) | 91.8 (90.2 to 93.4) | 92.9 (90.8 to 95.1) | 8 years | 9 years |
Fast Isokinetic Hamstrings
hs_fast_mv <- rma.mv(yi, hs_fastV,
mods = ~rcs(timepoint_mean, 4),
data = hs_fastdata,
random = list(~ timepoint_mean|cohort),
struct = c("CAR"))
robust(hs_fast_mv, cluster = cohort, clubSandwich = TRUE)##
## Multivariate Meta-Analysis Model (k = 79; method: REML)
##
## Variance Components:
##
## outer factor: cohort (nlvls = 53)
## inner factor: timepoint_mean (nlvls = 40)
##
## estim sqrt fixed
## tau^2 0.0094 0.0968 no
## rho 0.9105 no
##
## Test for Residual Heterogeneity:
## QE(df = 75) = 4054.9033, p-val < .0001
##
## Number of estimates: 79
## Number of clusters: 53
## Estimates per cluster: 1-6 (mean: 1.49, median: 1)
##
## Test of Moderators (coefficients 2:4):¹
## F(df1 = 3, df2 = 10.46) = 1.7558, p-val = 0.2163
##
## Model Results:
##
## estimate se¹ tval¹ df¹
## intrcpt -0.3249 0.1020 -3.1865 11.78
## rcs(timepoint_mean, 4)timepoint_mean 0.0401 0.0163 2.4573 12.87
## rcs(timepoint_mean, 4)timepoint_mean' -0.7001 0.2906 -2.4092 16.63
## rcs(timepoint_mean, 4)timepoint_mean'' 1.0612 0.4415 2.4036 17.04
## pval¹ ci.lb¹ ci.ub¹
## intrcpt 0.0080 -0.5476 -0.1023 **
## rcs(timepoint_mean, 4)timepoint_mean 0.0290 0.0048 0.0753 *
## rcs(timepoint_mean, 4)timepoint_mean' 0.0279 -1.3142 -0.0860 *
## rcs(timepoint_mean, 4)timepoint_mean'' 0.0279 0.1299 1.9925 *
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## 1) results based on cluster-robust inference (var-cov estimator: CR2,
## approx t/F-tests and confidence intervals, df: Satterthwaite approx)profile(hs_fast_mv)## Profiling tau2 = 1
## Profiling rho = 1
Plots - normal and log scale
mv_plotdetails(hs_fast_mv, include_pi = TRUE)
mv_plotdetails(hs_fast_mv, logscale = TRUE, showgraft = FALSE)
predict_details(hs_fast_mv) %>% kbl() %>% kable_styling(position = "left", full_width = FALSE)| 1 year | 2 years | 5 years | Zero crossing | Last Data Point |
|---|---|---|---|---|
| 94 (92 to 96) | 92.4 (90.1 to 94.7) | 95.6 (82.9 to 110.3) | 3.6 years | 5.1 years |
Isometric Hamstrings
hs_iso_mv <- rma.mv(yi, hs_isoV,
mods = ~log(timepoint_mean),
data = hs_isodata,
random = list(~ timepoint_mean|cohort),
struct = c("CAR"))
robust(hs_iso_mv, cluster = cohort, clubSandwich = TRUE)##
## Multivariate Meta-Analysis Model (k = 38; method: REML)
##
## Variance Components:
##
## outer factor: cohort (nlvls = 27)
## inner factor: timepoint_mean (nlvls = 31)
##
## estim sqrt fixed
## tau^2 0.0296 0.1720 no
## rho 0.7141 no
##
## Test for Residual Heterogeneity:
## QE(df = 36) = 4240.2591, p-val < .0001
##
## Number of estimates: 38
## Number of clusters: 27
## Estimates per cluster: 1-3 (mean: 1.41, median: 1)
##
## Test of Moderators (coefficient 2):¹
## F(df1 = 1, df2 = 14.82) = 1.4841, p-val = 0.2422
##
## Model Results:
##
## estimate se¹ tval¹ df¹ pval¹ ci.lb¹
## intrcpt -0.2836 0.1035 -2.7394 16.06 0.0145 -0.5030
## log(timepoint_mean) 0.0464 0.0381 1.2183 14.82 0.2422 -0.0349
## ci.ub¹
## intrcpt -0.0642 *
## log(timepoint_mean) 0.1278
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## 1) results based on cluster-robust inference (var-cov estimator: CR2,
## approx t/F-tests and confidence intervals, df: Satterthwaite approx)profile(hs_iso_mv)## Profiling tau2 = 1
## Profiling rho = 1
Plots - normal and log scale
mv_plotdetails(hs_iso_mv, include_pi = TRUE)
mv_plotdetails(hs_iso_mv, logscale = TRUE, showgraft = FALSE)
predict_details(hs_iso_mv) %>% kbl() %>% kable_styling(position = "left", full_width = FALSE)| 1 year | 2 years | 5 years | Zero crossing | Last Data Point |
|---|---|---|---|---|
| 84.5 (79.7 to 89.6) | 87.3 (80.9 to 94.2) | 91.1 (79.4 to 104.5) | 3.4 years | 5 years |
3. Quadriceps Strength - Case Control
# Generate the within person quad strength data
quadcont <- casecontrol %>%
filter(str_detect(measure, "quad")) %>%
mutate(es_id = row_number(),
timepoint_mean = acl_timepoint_mean,
group = 1) %>%
filter(timepoint_mean >2, # no pre-operative data
timepoint_mean < 120) %>%
mutate(measure_2 = case_when(
str_detect(measure, "mvic|hhd") ~ "Isometric",
str_detect(measure, "isk con 180|isk con 230|isk con 240|isk con 300") ~ "Fast Isokinetic",
str_detect(measure, "isk con 30|isk con 60|isk con 90|isk con 120") ~ "Slow Isokinetic"
)) %>%
filter(!is.na(measure_2)) %>%
#filter(timepoint_mean < 100) %>%
filter(!(cohort == "Laudner 2015" & measure == "quad isk con 300"), # remove where multiple effect sizes in same measure group
!(cohort == "Pamukoff 2018" & measure == "quad isk con 120"),
!(cohort == "Tourville 2014" & measure == "quad isk con 300"),
!(cohort == "Zult 2017" & measure == "quad isk con 120"))
# Reducing timepoint down to categories (3, 6, 12, 24, 48, 96 months post)
# not using this at the moment but keeping code here for now
quadcont_cat <- quadcont %>%
mutate(timepoint_cut = cut(timepoint_mean,
breaks = c(0, 4.5, 9, 18, 36, 72, Inf),
labels = c(3, 6, 12, 24, 48, 96)))
#group_by(cohort, timepoint_cut, group, measure_2) %>%
# if multiple timepoints allocated to same category, take the closest to the categorical timepoint
#slice(which.min(abs(timepoint_mean - as.numeric(as.character(timepoint_cut))))) %>%
#ungroup()
# separate subgroups of data based on contraction type
fastcont <- quadcont_cat %>% filter(measure_2 == "Fast Isokinetic")
slowcont <- quadcont_cat %>% filter(measure_2 == "Slow Isokinetic")
isocont <- quadcont_cat %>% filter(measure_2 == "Isometric")Specify covariance matrix
Need to create a full covariance matrix. Using clubSandwhich package to help Assuming a correlation of 0.85 between timepoints.
qcont_slowV <- impute_covariance_matrix(vi = slowcont$vi,
cluster = slowcont$cohort,
ti = slowcont$timepoint_mean,
ar1 = 0.85,
check_PD = TRUE,
return_list = FALSE)
qcont_fastV <- impute_covariance_matrix(vi = fastcont$vi,
cluster = fastcont$cohort,
ti = fastcont$timepoint_mean,
ar1 = 0.85,
check_PD = TRUE,
return_list = FALSE)
qcont_isoV <- impute_covariance_matrix(vi = isocont$vi,
cluster = isocont$cohort,
ti = isocont$timepoint_mean,
ar1 = 0.85,
check_PD = TRUE,
return_list = FALSE)Model Selection
Data are now ready for meta-analysis
# Create 'empty' models without timepoint used as a moderator variable.
# Will use these models to then build different versions on
# Using a continuous time auto-regressive variance structure (CAR)
qcont_slow_mv_empty <- rma.mv(yi, qcont_slowV,
data = slowcont,
random = list(~ timepoint_mean|cohort),
struct = c("CAR"))
qcont_fast_mv_empty <- rma.mv(yi, qcont_fastV,
data = fastcont,
random = list(~ timepoint_mean|cohort),
struct = c("CAR"))
qcont_iso_mv_empty <- rma.mv(yi, qcont_isoV,
data = isocont,
random = list(~ timepoint_mean|cohort),
struct = c("CAR"))Slow Isokinetic Quadriceps
l=Log best fit
mod_selection(qcont_slow_mv_empty)## [[1]]
## mod logLik. deviance. AIC. BIC. AICc.
## 1 Linear 12.06488 -24.12976 -16.12976 -11.765590 -13.776819
## 2 Log 13.34260 -26.68520 -18.68520 -14.321033 -16.332262
## 3 Poly (2) 11.23978 -22.47956 -12.47956 -7.256951 -8.479563
## 4 3 knot RCS 11.98789 -23.97579 -13.97579 -8.753174 -9.975787
## 5 4 knot RCS 13.01079 -26.02159 -14.02159 -8.047192 -7.560047
##
## [[2]]
Fast Isokinetic Quadriceps
Log spline best fit
mod_selection(qcont_fast_mv_empty)## [[1]]
## mod logLik. deviance. AIC. BIC. AICc.
## 1 Linear 8.183540 -16.36708 -8.3670806 -6.107283 -3.367081
## 2 Log 8.583040 -17.16608 -9.1660803 -6.906283 -4.166080
## 3 Poly (2) 7.366976 -14.73395 -4.7339514 -2.309418 5.266049
## 4 3 knot RCS 7.207663 -14.41533 -4.4153256 -1.990792 5.584674
## 5 4 knot RCS 6.281128 -12.56226 -0.5622562 1.825115 20.437744
##
## [[2]]
Isometric Quadriceps
Log is best fit
mod_selection(qcont_iso_mv_empty)## [[1]]
## mod logLik. deviance. AIC. BIC. AICc.
## 1 Linear 5.147059 -10.29412 -2.294118 2.8892299 -0.4759357
## 2 Log 5.988042 -11.97608 -3.976083 1.2072643 -2.1579014
## 3 Poly (2) 5.882500 -11.76500 -1.765000 4.5254823 1.2349997
## 4 3 knot RCS 5.602298 -11.20460 -1.204597 5.0858859 1.7954032
## 5 4 knot RCS 9.264707 -18.52941 -6.529415 0.7838403 -1.8627480
##
## [[2]]
Final models
Slow Isokinetic Quadriceps
qcont_slow_mv <- rma.mv(yi, qcont_slowV,
mods = ~log(timepoint_mean),
data = slowcont,
random = list(~ timepoint_mean|cohort),
struct = c("CAR"))
robust(qcont_slow_mv, cluster = cohort, clubSandwich = TRUE)##
## Multivariate Meta-Analysis Model (k = 24; method: REML)
##
## Variance Components:
##
## outer factor: cohort (nlvls = 18)
## inner factor: timepoint_mean (nlvls = 22)
##
## estim sqrt fixed
## tau^2 0.0152 0.1232 no
## rho 0.6874 no
##
## Test for Residual Heterogeneity:
## QE(df = 22) = 210.5326, p-val < .0001
##
## Number of estimates: 24
## Number of clusters: 18
## Estimates per cluster: 1-4 (mean: 1.33, median: 1)
##
## Test of Moderators (coefficient 2):¹
## F(df1 = 1, df2 = 7.87) = 17.6016, p-val = 0.0031
##
## Model Results:
##
## estimate se¹ tval¹ df¹ pval¹ ci.lb¹
## intrcpt -0.5194 0.0661 -7.8558 9.65 <.0001 -0.6675
## log(timepoint_mean) 0.1210 0.0288 4.1954 7.87 0.0031 0.0543
## ci.ub¹
## intrcpt -0.3714 ***
## log(timepoint_mean) 0.1877 **
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## 1) results based on cluster-robust inference (var-cov estimator: CR2,
## approx t/F-tests and confidence intervals, df: Satterthwaite approx)profile(qcont_slow_mv)## Profiling tau2 = 1
## Profiling rho = 1
Plots - normal and log scale
mv_plotdetails(qcont_slow_mv, include_pi = TRUE)
mv_plotdetails(qcont_slow_mv, logscale = TRUE, showgraft = FALSE)
predict_details(qcont_slow_mv) %>% kbl() %>% kable_styling(position = "left", full_width = FALSE)| 1 year | 2 years | 5 years | Zero crossing | Last Data Point |
|---|---|---|---|---|
| 80.3 (75.1 to 85.9) | 87.4 (78.8 to 96.9) | 97.6 (83.1 to 114.7) | 2.4 years | 4.6 years |
Fast Isokinetic Quadriceps
qcont_fast_mv <- rma.mv(yi, qcont_fastV,
mods = ~log(timepoint_mean),
data = fastcont,
random = list(~ timepoint_mean|cohort),
struct = c("CAR"))
robust(qcont_fast_mv, cluster = cohort, clubSandwich = TRUE)##
## Multivariate Meta-Analysis Model (k = 15; method: REML)
##
## Variance Components:
##
## outer factor: cohort (nlvls = 11)
## inner factor: timepoint_mean (nlvls = 13)
##
## estim sqrt fixed
## tau^2 0.0124 0.1112 no
## rho 0.0000 no
##
## Test for Residual Heterogeneity:
## QE(df = 13) = 116.5288, p-val < .0001
##
## Number of estimates: 15
## Number of clusters: 11
## Estimates per cluster: 1-4 (mean: 1.36, median: 1)
##
## Test of Moderators (coefficient 2):¹
## F(df1 = 1, df2 = 4.85) = 5.4999, p-val = 0.0675
##
## Model Results:
##
## estimate se¹ tval¹ df¹ pval¹ ci.lb¹
## intrcpt -0.3707 0.0664 -5.5795 4 0.0051 -0.5552
## log(timepoint_mean) 0.0580 0.0247 2.3452 4.85 0.0675 -0.0062
## ci.ub¹
## intrcpt -0.1862 **
## log(timepoint_mean) 0.1221 .
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## 1) results based on cluster-robust inference (var-cov estimator: CR2,
## approx t/F-tests and confidence intervals, df: Satterthwaite approx)profile(qcont_fast_mv)## Profiling tau2 = 1
## Profiling rho = 1
Plots - normal and log scale
mv_plotdetails(qcont_fast_mv, include_pi = TRUE)
mv_plotdetails(qcont_fast_mv, logscale = TRUE, showgraft = FALSE)
predict_details(qcont_fast_mv) %>% kbl() %>% kable_styling(position = "left", full_width = FALSE)| 1 year | 2 years | 5 years | Zero crossing | Last Data Point |
|---|---|---|---|---|
| 79.7 (73.6 to 86.3) | 83 (75.3 to 91.4) | 87.5 (75.8 to 101.1) | 4.5 years | 4.9 years |
Isometric Quadriceps
qcont_iso_mv <- rma.mv(yi, qcont_isoV,
mods = ~log(timepoint_mean),
data = isocont,
random = list(~ timepoint_mean|cohort),
struct = c("CAR"))
robust(qcont_iso_mv, cluster = cohort, clubSandwich = TRUE)##
## Multivariate Meta-Analysis Model (k = 29; method: REML)
##
## Variance Components:
##
## outer factor: cohort (nlvls = 27)
## inner factor: timepoint_mean (nlvls = 28)
##
## estim sqrt fixed
## tau^2 0.0337 0.1836 no
## rho 0.0000 no
##
## Test for Residual Heterogeneity:
## QE(df = 27) = 720.8710, p-val < .0001
##
## Number of estimates: 29
## Number of clusters: 27
## Estimates per cluster: 1-2 (mean: 1.07, median: 1)
##
## Test of Moderators (coefficient 2):¹
## F(df1 = 1, df2 = 12.03) = 5.1565, p-val = 0.0423
##
## Model Results:
##
## estimate se¹ tval¹ df¹ pval¹ ci.lb¹
## intrcpt -0.3777 0.0830 -4.5520 7.56 0.0022 -0.5710
## log(timepoint_mean) 0.0583 0.0257 2.2708 12.03 0.0423 0.0024
## ci.ub¹
## intrcpt -0.1844 **
## log(timepoint_mean) 0.1142 *
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## 1) results based on cluster-robust inference (var-cov estimator: CR2,
## approx t/F-tests and confidence intervals, df: Satterthwaite approx)profile(qcont_iso_mv)## Profiling tau2 = 1
## Profiling rho = 1
Plots - normal and log scale
mv_plotdetails(qcont_iso_mv, include_pi = TRUE)
mv_plotdetails(qcont_iso_mv, logscale = TRUE, showgraft = FALSE)
predict_details(qcont_iso_mv) %>% kbl() %>% kable_styling(position = "left", full_width = FALSE)| 1 year | 2 years | 5 years | Zero crossing | Last Data Point |
|---|---|---|---|---|
| 79.2 (73.9 to 85) | 82.5 (77.6 to 87.7) | 87 (79.9 to 94.7) | 8.7 years | 4.5 years |
4. Hamstring Strength - Case Control
hscont <- casecontrol %>%
filter(str_detect(measure, "hs")) %>%
mutate(es_id = row_number(),
timepoint_mean = acl_timepoint_mean,
group = 1) %>%
filter(timepoint_mean >2, # no pre-operative data
timepoint_mean < 120) %>%
mutate(measure_2 = case_when(
str_detect(measure, "mvic|hhd") ~ "Isometric",
str_detect(measure, "isk con 180|isk con 230|isk con 240|isk con 300") ~ "Fast Isokinetic",
str_detect(measure, "isk con 30|isk con 60|isk con 90|isk con 120") ~ "Slow Isokinetic"
)) %>%
filter(!is.na(measure_2)) %>%
filter(timepoint_mean < 100) %>%
filter(!(cohort == "Laudner 2015" & measure == "hs isk con 300"), # remove where multiple effect sizes in same measure group
!(cohort == "Tourville 2014" & measure == "hs isk con 300"),
!(cohort == "Zult 2017" & measure == "hs isk con 120"))
# Reducing timepoint down to categories (3, 6, 12, 24, 48, 96 months post)
# not using this at the moment but keeping code here for now
hscont_cat <- hscont %>%
mutate(timepoint_cut = cut(timepoint_mean,
breaks = c(0, 4.5, 9, 18, 36, 72, Inf),
labels = c(3, 6, 12, 24, 48, 96)))
#group_by(cohort, timepoint_cut, measure_2) %>%
# if multiple timepoints allocated to same category, take the closest to the categorical timepoint
#slice(which.min(abs(timepoint_mean - as.numeric(as.character(timepoint_cut))))) %>%
#ungroup()
# separate subgroups of data based on contraction type
hs_fastcont <- hscont_cat %>% filter(measure_2 == "Fast Isokinetic")
hs_slowcont <- hscont_cat %>% filter(measure_2 == "Slow Isokinetic")
hs_isocont <- hscont_cat %>% filter(measure_2 == "Isometric")Specify covariance matrix
Need to create a full covariance matrix. Using clubSandwhich package to help Assuming a correlation of 0.85 between timepoints.
hcont_fastV <- impute_covariance_matrix(vi = hs_fastcont$vi,
cluster = hs_fastcont$cohort,
ti = hs_fastcont$timepoint_mean,
ar1 = 0.85,
check_PD = TRUE,
smooth_vi = TRUE,
return_list = FALSE)
hcont_slowV <- impute_covariance_matrix(vi = hs_slowcont$vi,
cluster = hs_slowcont$cohort,
ti = hs_slowcont$timepoint_mean,
ar1 = 0.85,
smooth_vi = TRUE,
return_list = FALSE)
hcont_isoV <- impute_covariance_matrix(vi = hs_isocont$vi,
cluster = hs_isocont$cohort,
ti = hs_isocont$timepoint_mean,
ar1 = 0.85,
check_PD = TRUE,
smooth_vi = TRUE,
return_list = FALSE)Model Selection
# Create 'empty' models without timepoint used as a moderator variable.
# Will use these models to then build different versions on
# Using a continuous time auto-regressive variance structure (CAR)
hcont_slow_mv_empty <- rma.mv(yi, hcont_slowV,
data = hs_slowcont,
random = list(~ timepoint_mean|cohort),
struct = c("CAR"))
hcont_fast_mv_empty <- rma.mv(yi, hcont_fastV,
data = hs_fastcont,
random = list(~ timepoint_mean|cohort),
struct = c("CAR"))
hcont_iso_mv_empty <- rma.mv(yi, hcont_isoV,
data = hs_isocont,
random = list(~ timepoint_mean|cohort),
struct = c("CAR"))Slow Isokinetic Hamstrings
Log best fit
mod_selection(hcont_slow_mv_empty)## [[1]]
## mod logLik. deviance. AIC. BIC. AICc.
## 1 Linear 13.20571 -26.41141 -18.41141 -15.32106 -14.775050
## 2 Log 15.04066 -30.08131 -22.08131 -18.99096 -18.444947
## 3 Poly (2) 12.31979 -24.63957 -14.63957 -11.09932 -7.972908
## 4 3 knot RCS 16.94020 -33.88040 -23.88040 -20.34015 -17.213733
## 5 4 knot RCS 18.36443 -36.72885 -24.72885 -20.89451 -12.728854
##
## [[2]]
Fast Isokinetic Hamstrings
Log best fit
mod_selection(hcont_fast_mv_empty)## [[1]]
## mod logLik. deviance. AIC. BIC. AICc.
## 1 Linear 3.937956 -7.875911 0.1240888 1.715670 6.790755
## 2 Log 4.052934 -8.105868 -0.1058677 1.485713 6.560799
## 3 Poly (2) 3.162075 -6.324150 3.6758500 5.188775 18.675850
## 4 3 knot RCS 3.483855 -6.967710 3.0322903 4.545216 18.032290
## 5 4 knot RCS 3.827806 -7.655612 2.3443884 3.857314 17.344388
##
## [[2]]
Isometric hamstrings
Polynomial (2) is best fit A lot of uncertainty here.
mod_selection(hcont_iso_mv_empty)## [[1]]
## mod logLik. deviance. AIC. BIC. AICc.
## 1 Linear 2.465917 -4.931833 3.0681665 5.007793 8.782452
## 2 Log 1.821243 -3.642486 4.3575138 6.297140 10.071800
## 3 Poly (2) 4.590334 -9.180668 0.8193318 2.808808 12.819332
## 4 3 knot RCS 4.417871 -8.835742 1.1642579 3.153734 13.164258
## 5 4 knot RCS 5.223635 -10.447270 1.5527296 3.368240 29.552730
##
## [[2]]
Final models
Slow Isokinetic Hamstrings
hcont_slow_mv <- rma.mv(yi, hcont_slowV,
mods = ~log(timepoint_mean),
data = hs_slowcont,
random = list(~ timepoint_mean|cohort),
struct = c("CAR"))
robust(hcont_slow_mv, cluster = cohort, clubSandwich = TRUE)##
## Multivariate Meta-Analysis Model (k = 18; method: REML)
##
## Variance Components:
##
## outer factor: cohort (nlvls = 12)
## inner factor: timepoint_mean (nlvls = 17)
##
## estim sqrt fixed
## tau^2 0.0092 0.0962 no
## rho 0.8840 no
##
## Test for Residual Heterogeneity:
## QE(df = 16) = 97.2809, p-val < .0001
##
## Number of estimates: 18
## Number of clusters: 12
## Estimates per cluster: 1-4 (mean: 1.50, median: 1)
##
## Test of Moderators (coefficient 2):¹
## F(df1 = 1, df2 = 4.32) = 23.6758, p-val = 0.0068
##
## Model Results:
##
## estimate se¹ tval¹ df¹ pval¹ ci.lb¹
## intrcpt -0.2900 0.0624 -4.6463 6.36 0.0030 -0.4407
## log(timepoint_mean) 0.0983 0.0202 4.8658 4.32 0.0068 0.0438
## ci.ub¹
## intrcpt -0.1393 **
## log(timepoint_mean) 0.1528 **
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## 1) results based on cluster-robust inference (var-cov estimator: CR2,
## approx t/F-tests and confidence intervals, df: Satterthwaite approx)profile(hcont_slow_mv)## Profiling tau2 = 1
## Profiling rho = 1
Plots - normal and log scale
mv_plotdetails(hcont_slow_mv, include_pi = TRUE)
mv_plotdetails(hcont_slow_mv, logscale = TRUE, showgraft = FALSE)
predict_details(hcont_slow_mv) %>% kbl() %>% kable_styling(position = "left", full_width = FALSE)| 1 year | 2 years | 5 years | Zero crossing | Last Data Point |
|---|---|---|---|---|
| 95.5 (88.7 to 102.9) | 102.3 (92.7 to 112.8) | 111.9 (97.7 to 128.1) | 0.7 years | 4.6 years |
Fast Isokinetic Hamstrings
hcont_fast_mv <- rma.mv(yi, hcont_fastV,
mods = ~log(timepoint_mean),
data = hs_fastcont,
random = list(~ timepoint_mean|cohort),
struct = c("CAR"))
robust(hcont_fast_mv, cluster = cohort, clubSandwich = TRUE)##
## Multivariate Meta-Analysis Model (k = 13; method: REML)
##
## Variance Components:
##
## outer factor: cohort (nlvls = 9)
## inner factor: timepoint_mean (nlvls = 12)
##
## estim sqrt fixed
## tau^2 0.0326 0.1804 no
## rho 0.9154 no
##
## Test for Residual Heterogeneity:
## QE(df = 11) = 102.0693, p-val < .0001
##
## Number of estimates: 13
## Number of clusters: 9
## Estimates per cluster: 1-4 (mean: 1.44, median: 1)
##
## Test of Moderators (coefficient 2):¹
## F(df1 = 1, df2 = 2.56) = 19.9575, p-val = 0.0290
##
## Model Results:
##
## estimate se¹ tval¹ df¹ pval¹ ci.lb¹
## intrcpt -0.3128 0.0693 -4.5148 3.29 0.0166 -0.5227
## log(timepoint_mean) 0.0900 0.0202 4.4674 2.56 0.0290 0.0192
## ci.ub¹
## intrcpt -0.1030 *
## log(timepoint_mean) 0.1609 *
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## 1) results based on cluster-robust inference (var-cov estimator: CR2,
## approx t/F-tests and confidence intervals, df: Satterthwaite approx)profile(hcont_fast_mv)## Profiling tau2 = 1
## Profiling rho = 1
Plots - normal and log scale
mv_plotdetails(hcont_fast_mv, include_pi = TRUE)
mv_plotdetails(hcont_fast_mv, logscale = TRUE, showgraft = FALSE)
predict_details(hcont_fast_mv) %>% kbl() %>% kable_styling(position = "left", full_width = FALSE)| 1 year | 2 years | 5 years | Zero crossing | Last Data Point |
|---|---|---|---|---|
| 91.5 (79.5 to 105.3) | 97.4 (79.3 to 119.6) | 105.7 (80.3 to 139.3) | 0.5 years | 4.3 years |
Isometric Hamstrings
hcont_iso_mv <- rma.mv(yi, hcont_isoV,
mods = ~poly(timepoint_mean, degree = 2, raw = TRUE),
data = hs_isocont,
random = list(~ timepoint_mean|cohort),
struct = c("CAR"))
robust(hcont_iso_mv, cluster = cohort, clubSandwich = TRUE)##
## Multivariate Meta-Analysis Model (k = 14; method: REML)
##
## Variance Components:
##
## outer factor: cohort (nlvls = 12)
## inner factor: timepoint_mean (nlvls = 13)
##
## estim sqrt fixed
## tau^2 0.0238 0.1544 no
## rho 0.6188 no
##
## Test for Residual Heterogeneity:
## QE(df = 11) = 168.2278, p-val < .0001
##
## Number of estimates: 14
## Number of clusters: 12
## Estimates per cluster: 1-2 (mean: 1.17, median: 1)
##
## Test of Moderators (coefficients 2:3):¹
## F(df1 = 2, df2 = 4.37) = 15.8317, p-val = 0.0100
##
## Model Results:
##
## estimate se¹ tval¹
## intrcpt -0.3123 0.1578 -1.9789
## poly(timepoint_mean, degree = 2, raw = TRUE)1 0.0433 0.0170 2.5415
## poly(timepoint_mean, degree = 2, raw = TRUE)2 -0.0013 0.0004 -3.3905
## df¹ pval¹ ci.lb¹
## intrcpt 4.81 0.1070 -0.7229
## poly(timepoint_mean, degree = 2, raw = TRUE)1 6.05 0.0437 0.0017
## poly(timepoint_mean, degree = 2, raw = TRUE)2 5.61 0.0163 -0.0022
## ci.ub¹
## intrcpt 0.0984
## poly(timepoint_mean, degree = 2, raw = TRUE)1 0.0848 *
## poly(timepoint_mean, degree = 2, raw = TRUE)2 -0.0003 *
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## 1) results based on cluster-robust inference (var-cov estimator: CR2,
## approx t/F-tests and confidence intervals, df: Satterthwaite approx)profile(hcont_iso_mv)## Profiling tau2 = 1
## Profiling rho = 1
Plots - normal and log scale
mv_plotdetails(hcont_iso_mv, include_pi = TRUE)
mv_plotdetails(hcont_iso_mv, logscale = TRUE, showgraft = FALSE)
predict_details(hcont_iso_mv) %>% kbl() %>% kable_styling(position = "left", full_width = FALSE)| 1 year | 2 years | 5 years | Zero crossing | Last Data Point |
|---|---|---|---|---|
| 102.1 (89.3 to 116.7) | 98.1 (84.8 to 113.5) | 9.3 (2.4 to 35.5) | 2.3 years | 3 years |
5. Bivariate Analysis - Comparing Between and Within person comparisons
Of interest is whether studies that compare both between person and within person show similar effects - i.e. are the effects different For this we select all the datapoints that have both between and within person comparisons and conduct a bivariate meta-analysis. We have to account for non-independence of samples in the analysis conducting a variance/covariance matrix and also fitting random effects for this (type of outcome: within/between person nested within cohorts). To calculate this we assumed a rho of 0.7. Sensitivity analyses show that estimates are stable to different values for this.
Because the majority of studies report only between person comparisons, we chose to only look at those studies that measured both (i.e. “bivariate”), otherwise estimation of the correlation between effects will be spurious.
Set up:
## Bivariate Analyses
## Investigating relationship between within and between person estimates
# Within Person Data
quad1 <- quad %>%
select(study, cohort, measure, measure_2, acl_n, timepoint, timepoint_mean, group, acl_graft_group, yi, vi) %>%
rename(acl_group = group) %>%
mutate(type = "within")
# Between person data
quad2 <- quadcont %>%
select(study, cohort, measure, measure_2, acl_n, timepoint, timepoint_mean, acl_group, acl_graft_group, yi, vi) %>%
mutate(type = "casecontrol")
# Join together
quad_joined <- bind_rows(quad1, quad2)
# Split based on contraction types
slow_joined <- quad_joined %>% filter(measure_2 == "Slow Isokinetic") %>%
filter(!(cohort == "Li 2022" & measure == "quad isk con 30"),
!(cohort == "Siney 2010" & measure == "quad isk con 90")) %>%
group_by(study, measure, timepoint) %>%
filter(length(type) == 2) %>% # make sure to only include studies who measure both effects
ungroup() %>%
group_by(study, measure, type) %>%
slice(which.min(abs(timepoint_mean - as.numeric(12)))) %>% # take one effect per study, preferably closest to 12 months if multiple
ungroup() %>%
filter(!is.na(vi)) %>% # remove any missing data
group_by(study) %>%
filter(n() > 1) %>%
ungroup()
fast_joined <- quad_joined %>% filter(measure_2 == "Fast Isokinetic") %>%
filter(!(cohort == "Bailey 2019" & measure == "quad isk con 300"),
!(cohort == "Drocco 2017" & measure == "quad isk con 240"),
!(cohort == "Kyritsis 2016" & measure == "quad isk con 300"),
!(cohort == "Laudner 2015" & measure == "quad isk con 300"),
!(cohort == "Pamukoff 2018" & measure == "quad isk con 120"),
!(cohort == "Tourville 2014" & measure == "quad isk con 300"),
!(cohort == "Welling 2020" & measure == "quad isk con 300"),
!(cohort == "Welling 2018b" & measure == "quad isk con 300"),
!(cohort == "Zult 2017" & measure == "quad isk con 120")) %>%
group_by(study, measure, timepoint) %>%
filter(length(type) == 2) %>%
ungroup() %>%
group_by(study, measure, type) %>%
slice(which.min(abs(timepoint_mean - as.numeric(12)))) %>%
ungroup() %>%
filter(!is.na(vi)) %>% # remove any missing data
group_by(study) %>%
filter(n() > 1) %>%
ungroup()
iso_joined <- quad_joined %>% filter(measure_2 == "Isometric") %>%
filter(!(cohort == "Karanikas 2005" & measure == "quad mvic 0"),
!(cohort == "Labanca 2018" & measure == "quad mvic 30"),
!(cohort == "Labanca 2016" & measure == "quad mvic 30")) %>%
group_by(study, measure, timepoint) %>%
filter(length(type) == 2) %>%
ungroup() %>%
group_by(study, measure, type) %>%
slice(which.min(abs(timepoint_mean - as.numeric(12)))) %>%
ungroup() %>%
filter(!is.na(vi)) %>% # remove any missing data
group_by(study) %>%
filter(n() > 1) %>%
ungroup()
# Create covariance matrices for each, assuming rho of 0.7
slowjoinedV <- vcalc(vi,
cluster = study,
type = type,
rho = 0.7,
data = slow_joined,
checkpd = TRUE)
fastjoinedV <- vcalc(vi,
cluster = study,
type = type,
rho = 0.7,
data = fast_joined,
checkpd = TRUE)
isojoinedV <- vcalc(vi,
cluster = study,
type = type,
rho = 0.7,
data = iso_joined,
checkpd = TRUE)
# Fit individual models with random effects for type of effect, nested within each cohort.
slow_joined_mv <- rma.mv(yi,
slowjoinedV,
mods = ~ type - 1, # remove intercept to generate estimate for each.
random = list(~ type | cohort),
struct = c("UN"),
data = slow_joined,
cvvc = "varcov") # need this to use later with matreg
fast_joined_mv <- rma.mv(yi,
fastjoinedV,
mods = ~ type - 1,
random = list(~ type | cohort),
struct = c("UN"),
data = fast_joined,
cvvc = "varcov")
iso_joined_mv <- rma.mv(yi,
isojoinedV,
mods = ~ type - 1,
random = list(~ type | cohort),
struct = c("UN"),
data = iso_joined,
cvvc = "varcov")
####
## Hamstring
####
hs1 <- hs %>%
select(study, cohort, measure, measure_2, acl_n, timepoint, timepoint_mean, group, acl_graft_group, yi, vi) %>%
rename(acl_group = group) %>%
mutate(type = "within")
hs2 <- hscont %>%
select(study, cohort, measure, measure_2, acl_n, timepoint, timepoint_mean, acl_group, acl_graft_group, yi, vi) %>%
mutate(type = "casecontrol")
hs_joined <- bind_rows(hs1, hs2)
slow_joined_hs <- hs_joined %>% filter(measure_2 == "Slow Isokinetic") %>%
filter(!(cohort == "Lee 2019" & measure == "hs isk con 60 - deep"),
!(cohort == "Kim 2011" & measure == "hs isk con 90 hyperflex"),
!(cohort == "Siney 2010" & measure == "hs isk con 90"),
!(cohort == "McRae 2013" & measure == "hs isk con 60 supine"),
!(cohort == "Li 2022" & measure == "hs isk con 30")) %>%
group_by(study, measure, timepoint) %>%
filter(length(type) == 2) %>%
ungroup() %>%
group_by(study, measure, type) %>%
slice(which.min(abs(timepoint_mean - as.numeric(12)))) %>%
ungroup() %>%
filter(!is.na(vi)) %>% # remove any missing data
group_by(study) %>%
filter(n() > 1) %>%
ungroup()
fast_joined_hs <- hs_joined %>% filter(measure_2 == "Fast Isokinetic") %>%
filter(!(cohort == "Drocco 2017" & measure == "hs isk con 240"),
!(cohort == "Kyritsis 2016" & measure == "hs isk con 300"),
!(cohort == "Laudner 2015" & measure == "hs isk con 300"),
!(cohort == "Tourville 2014" & measure == "hs isk con 300"),
!(cohort == "Welling 2020" & measure == "hs isk con 300"),
!(cohort == "Welling 2018b" & measure == "hs isk con 300")) %>%
group_by(study, measure, timepoint) %>%
filter(length(type) == 2) %>%
ungroup() %>%
group_by(study, measure, type) %>%
slice(which.min(abs(timepoint_mean - as.numeric(12)))) %>%
ungroup() %>%
filter(!is.na(vi)) %>% # remove any missing data
group_by(study) %>%
filter(n() > 1) %>%
ungroup()
iso_joined_hs <- hs_joined %>% filter(measure_2 == "Isometric") %>%
filter(!(cohort == "Hu 2020" & measure == "hs mvic 30 hhd"),
!(cohort == "Hu 2020" & measure == "hs mvic 60 hhd"),
!(cohort == "Ardern 2010" & measure == "hs mvic 30"),
!(cohort == "Ardern 2010" & measure == "hs mvic 105")) %>%
group_by(study, measure, timepoint) %>%
filter(length(type) == 2) %>%
ungroup() %>%
group_by(study, measure, type) %>%
slice(which.min(abs(timepoint_mean - as.numeric(12)))) %>%
ungroup() %>%
filter(!is.na(vi)) %>% # remove any missing data
group_by(study) %>%
filter(n() > 1) %>%
ungroup()
slowjoined_hsV <- vcalc(vi,
cluster = study,
type = type,
rho = 0.7,
data = slow_joined_hs,
checkpd = TRUE)
fastjoined_hsV <- vcalc(vi,
cluster = study,
type = type,
rho = 0.7,
data = fast_joined_hs,
checkpd = TRUE)
isojoined_hsV <- vcalc(vi,
cluster = study,
type = type,
rho = 0.7,
data = iso_joined_hs,
checkpd = TRUE)
slow_joined_hs_mv <- rma.mv(yi,
slowjoined_hsV,
mods = ~ type - 1,
random = list(~ type | cohort),
struct = c("UN"),
data = slow_joined_hs,
cvvc = "varcov")
fast_joined_hs_mv <- rma.mv(yi,
fastjoined_hsV,
mods = ~ type - 1,
random = list(~ type | cohort),
struct = c("UN"),
data = fast_joined_hs,
cvvc = "varcov")
iso_joined_hs_mv <- rma.mv(yi,
isojoined_hsV,
mods = ~ type - 1,
random = list(~ type | cohort),
struct = c("UN"),
data = iso_joined_hs,
cvvc = "varcov")For model simplicity we only selected one estimate per study (i.e. one timepoint), as almost all between person studies only measure one timepoint. For the few studies that report multiple, we took the timepoint closest to 12 months (most common timepoint)
The resulting analysis provides an estimate for both between and within person estimates. More interestingly we also get an estimate of the correlation between effects.
We can then regress the true effects from the model to see the relationship between the two outcomes
In general we can say that between person effects don’t always mirror within person effects
Quadriceps - Slow Isokinetic
summary(slow_joined_mv)##
## Multivariate Meta-Analysis Model (k = 32; method: REML)
##
## logLik Deviance AIC BIC AICc
## 18.2483 -36.4967 -26.4967 -19.4907 -23.9967
##
## Variance Components:
##
## outer factor: cohort (nlvls = 16)
## inner factor: type (nlvls = 2)
##
## estim sqrt k.lvl fixed level
## tau^2.1 0.0277 0.1665 16 no casecontrol
## tau^2.2 0.0181 0.1347 16 no within
##
## rho.cscn rho.wthn cscn wthn
## casecontrol 1 - 16
## within 0.7132 1 no -
##
## Test for Residual Heterogeneity:
## QE(df = 30) = 969.4979, p-val < .0001
##
## Test of Moderators (coefficients 1:2):
## QM(df = 2) = 30.9018, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## typecasecontrol -0.2124 0.0433 -4.9016 <.0001 -0.2973 -0.1275 ***
## typewithin -0.1811 0.0341 -5.3171 <.0001 -0.2478 -0.1143 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1corplot_function(slow_joined_mv)## [[1]]
##
## estimate se zval pval ci.lb ci.ub
## intrcpt 0.3681 0.1312 2.8055 0.0050 0.1109 0.6252 **
## casecontrol 0.5767 0.1622 3.5546 0.0004 0.2587 0.8946 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## [[2]]
Quadriceps - Fast Isokinetic
summary(fast_joined_mv)##
## Multivariate Meta-Analysis Model (k = 20; method: REML)
##
## logLik Deviance AIC BIC AICc
## 15.0261 -30.0522 -20.0522 -15.6003 -15.0522
##
## Variance Components:
##
## outer factor: cohort (nlvls = 10)
## inner factor: type (nlvls = 2)
##
## estim sqrt k.lvl fixed level
## tau^2.1 0.0125 0.1119 10 no casecontrol
## tau^2.2 0.0080 0.0894 10 no within
##
## rho.cscn rho.wthn cscn wthn
## casecontrol 1 - 10
## within 0.4530 1 no -
##
## Test for Residual Heterogeneity:
## QE(df = 18) = 298.7612, p-val < .0001
##
## Test of Moderators (coefficients 1:2):
## QM(df = 2) = 34.7953, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## typecasecontrol -0.1969 0.0396 -4.9678 <.0001 -0.2746 -0.1192 ***
## typewithin -0.1499 0.0292 -5.1350 <.0001 -0.2071 -0.0927 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1corplot_function(fast_joined_mv)## [[1]]
##
## estimate se zval pval ci.lb ci.ub
## intrcpt 0.5636 0.2420 2.3287 0.0199 0.0892 1.0379 *
## casecontrol 0.3619 0.2947 1.2282 0.2194 -0.2156 0.9395
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## [[2]]
Quadriceps - Isometric
summary(iso_joined_mv)##
## Multivariate Meta-Analysis Model (k = 28; method: REML)
##
## logLik Deviance AIC BIC AICc
## 14.8694 -29.7387 -19.7387 -13.4483 -16.7387
##
## Variance Components:
##
## outer factor: cohort (nlvls = 14)
## inner factor: type (nlvls = 2)
##
## estim sqrt k.lvl fixed level
## tau^2.1 0.0340 0.1843 14 no casecontrol
## tau^2.2 0.0238 0.1542 14 no within
##
## rho.cscn rho.wthn cscn wthn
## casecontrol 1 - 14
## within 0.7957 1 no -
##
## Test for Residual Heterogeneity:
## QE(df = 26) = 830.8786, p-val < .0001
##
## Test of Moderators (coefficients 1:2):
## QM(df = 2) = 18.5663, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## typecasecontrol -0.2103 0.0517 -4.0707 <.0001 -0.3116 -0.1091 ***
## typewithin -0.1692 0.0417 -4.0597 <.0001 -0.2508 -0.0875 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1corplot_function(iso_joined_mv)## [[1]]
##
## estimate se zval pval ci.lb ci.ub
## intrcpt 0.3051 0.1234 2.4719 0.0134 0.0632 0.5470 *
## casecontrol 0.6655 0.1523 4.3692 <.0001 0.3670 0.9640 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## [[2]]
Hamstrings - Slow Isokinetic
summary(slow_joined_hs_mv)##
## Multivariate Meta-Analysis Model (k = 20; method: REML)
##
## logLik Deviance AIC BIC AICc
## 27.8370 -55.6741 -45.6741 -41.2222 -40.6741
##
## Variance Components:
##
## outer factor: cohort (nlvls = 10)
## inner factor: type (nlvls = 2)
##
## estim sqrt k.lvl fixed level
## tau^2.1 0.0062 0.0787 10 no casecontrol
## tau^2.2 0.0007 0.0269 10 no within
##
## rho.cscn rho.wthn cscn wthn
## casecontrol 1 - 10
## within 0.2382 1 no -
##
## Test for Residual Heterogeneity:
## QE(df = 18) = 128.1188, p-val < .0001
##
## Test of Moderators (coefficients 1:2):
## QM(df = 2) = 32.7972, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## typecasecontrol -0.0215 0.0280 -0.7676 0.4427 -0.0763 0.0333
## typewithin -0.0548 0.0098 -5.5996 <.0001 -0.0740 -0.0356 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1corplot_function(slow_joined_hs_mv)## [[1]]
##
## estimate se zval pval ci.lb ci.ub
## intrcpt 0.8669 0.1422 6.0956 <.0001 0.5882 1.1457 ***
## casecontrol 0.0815 0.1453 0.5607 0.5750 -0.2033 0.3663
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## [[2]]
Hamstrings - Fast Isokinetic
summary(fast_joined_hs_mv)##
## Multivariate Meta-Analysis Model (k = 14; method: REML)
##
## logLik Deviance AIC BIC AICc
## 11.0240 -22.0480 -12.0480 -9.6235 -2.0480
##
## Variance Components:
##
## outer factor: cohort (nlvls = 7)
## inner factor: type (nlvls = 2)
##
## estim sqrt k.lvl fixed level
## tau^2.1 0.0459 0.2144 7 no casecontrol
## tau^2.2 0.0061 0.0780 7 no within
##
## rho.cscn rho.wthn cscn wthn
## casecontrol 1 - 7
## within 0.8629 1 no -
##
## Test for Residual Heterogeneity:
## QE(df = 12) = 167.9627, p-val < .0001
##
## Test of Moderators (coefficients 1:2):
## QM(df = 2) = 6.4110, p-val = 0.0405
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## typecasecontrol -0.0686 0.0832 -0.8239 0.4100 -0.2316 0.0945
## typewithin -0.0593 0.0304 -1.9506 0.0511 -0.1188 0.0003 .
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1corplot_function(fast_joined_hs_mv)## [[1]]
##
## estimate se zval pval ci.lb ci.ub
## intrcpt 0.6493 0.0850 7.6421 <.0001 0.4827 0.8158 ***
## casecontrol 0.3140 0.0910 3.4512 0.0006 0.1357 0.4923 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## [[2]]
Hamstrings - Isometric
summary(iso_joined_hs_mv)##
## Multivariate Meta-Analysis Model (k = 16; method: REML)
##
## logLik Deviance AIC BIC AICc
## 13.6682 -27.3365 -17.3365 -14.1412 -9.8365
##
## Variance Components:
##
## outer factor: cohort (nlvls = 8)
## inner factor: type (nlvls = 2)
##
## estim sqrt k.lvl fixed level
## tau^2.1 0.0353 0.1879 8 no casecontrol
## tau^2.2 0.0016 0.0401 8 no within
##
## rho.cscn rho.wthn cscn wthn
## casecontrol 1 - 8
## within 0.5220 1 no -
##
## Test for Residual Heterogeneity:
## QE(df = 14) = 170.0779, p-val < .0001
##
## Test of Moderators (coefficients 1:2):
## QM(df = 2) = 39.2397, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## typecasecontrol -0.0902 0.0698 -1.2923 0.1962 -0.2270 0.0466
## typewithin -0.0965 0.0164 -5.8691 <.0001 -0.1288 -0.0643 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1corplot_function(iso_joined_hs_mv)## [[1]]
##
## estimate se zval pval ci.lb ci.ub
## intrcpt 0.8063 0.0766 10.5209 <.0001 0.6561 0.9565 ***
## casecontrol 0.1113 0.0839 1.3273 0.1844 -0.0531 0.2757
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## [[2]]
6. Summary
Quadriceps - Slow
Quadriceps - Fast
Quadriceps - Isometric
Hamstrings - Slow
Hamstrings - Fast
Hamstrings - Isometric
7. Eccentric Contractinos
ecc_within <- within_data %>%
filter(!study %in% c("Tengman 2014b", "Hogberg 2023")) %>%
rename(acl_timepoint_mean = timepoint_mean) %>%
mutate(n_acl_2 = as.character(as.integer(acl_n))) %>%
filter(str_detect(measure, "ecc"),
str_detect(measure, "quad|hs")) %>%
mutate(measure_2 = case_when(
str_detect(measure, "150|180|230|240|300") ~ "Fast Isokinetic",
str_detect(measure, "30|60|90") ~ "Slow Isokinetic"
)) %>%
mutate(measure_3 = case_when(
acl_timepoint_mean <= 12 & measure_2 == "Slow Isokinetic" ~ "Slow Isokinetic <= 1 year",
acl_timepoint_mean > 12 & measure_2 == "Slow Isokinetic" ~ "Slow Isokinetic > 1 year",
TRUE ~ measure_2),
measure_3 = factor(measure_3, levels = c("Slow Isokinetic <= 1 year", "Slow Isokinetic > 1 year",
"Fast Isokinetic"))
)
ecc_within$sei[ecc_within$study == "Gauthier 2022"] <- 0.01209894 # missing for variance for one study, taking value frmo other similar study Raoul.Quadriceps Eccentric - Within Person
Only 6 studies report eccentric quadriceps outcomes, most for slow eccentric, at variety of timepoints.
library(meta) # use meta package to fit univariate meta-analyses (for nicer plots)
quad_ecc_meta <- ecc_within %>%
filter(str_detect(measure, "quad")) %>%
metagen(TE = yi, seTE = sei, studlab = study, data = ., subgroup = measure_2, sm = "ROM")
summary(quad_ecc_meta)## ROM 95%-CI %W(common) %W(random)
## Harput 2018 0.8287 [0.7905; 0.8687] 19.9 24.7
## Engelen-VanMelick 2017 0.9844 [0.9521; 1.0179] 39.4 25.9
## Heijne 2015 0.9300 [0.8901; 0.9716] 23.0 25.0
## Heijne 2015 0.9600 [0.9132; 1.0092] 17.7 24.4
## measure_2
## Harput 2018 Slow Isokinetic
## Engelen-VanMelick 2017 Slow Isokinetic
## Heijne 2015 Slow Isokinetic
## Heijne 2015 Fast Isokinetic
##
## Number of studies combined: k = 4
##
## ROM 95%-CI z p-value
## Common effect model 0.9348 [0.9154; 0.9546] -6.29 < 0.0001
## Random effects model 0.9245 [0.8582; 0.9959] -2.07 0.0386
##
## Quantifying heterogeneity:
## tau^2 = 0.0053 [0.0014; 0.0800]; tau = 0.0726 [0.0369; 0.2829]
## I^2 = 91.5% [81.5%; 96.1%]; H = 3.44 [2.32; 5.08]
##
## Test of heterogeneity:
## Q d.f. p-value
## 35.44 3 < 0.0001
##
## Results for subgroups (common effect model):
## k ROM 95%-CI Q I^2
## measure_2 = Slow Isokinetic 3 0.9295 [0.9082; 0.9512] 34.12 94.1%
## measure_2 = Fast Isokinetic 1 0.9600 [0.9132; 1.0092] 0.00 --
##
## Test for subgroup differences (common effect model):
## Q d.f. p-value
## Between groups 1.33 1 0.2496
## Within groups 34.12 2 < 0.0001
##
## Results for subgroups (random effects model):
## k ROM 95%-CI tau^2 tau
## measure_2 = Slow Isokinetic 3 0.9130 [0.8268; 1.0082] 0.0072 0.0850
## measure_2 = Fast Isokinetic 1 0.9600 [0.9132; 1.0092] -- --
##
## Test for subgroup differences (random effects model):
## Q d.f. p-value
## Between groups 0.79 1 0.3754
##
## Details on meta-analytical method:
## - Inverse variance method
## - Restricted maximum-likelihood estimator for tau^2
## - Q-Profile method for confidence interval of tau^2 and tauforest_cc_function(quad_ecc_meta, overall = FALSE, prediction.subgroup = TRUE, weight.study = "same", hetstat = FALSE)
Hamstring Eccentric - Within Person
11 studies report eccentric quadriceps outcomes, most for slow eccentric, at variety of timepoints.
hs_ecc_meta <- ecc_within %>%
filter(str_detect(measure, "hs")) %>%
metagen(TE = yi, seTE = sei, studlab = study, data = ., subgroup = measure_3, sm = "ROM")
summary(hs_ecc_meta)## ROM 95%-CI %W(common) %W(random)
## Joseph 2020 0.8848 [0.8539; 0.9168] 5.9 9.4
## Raoul 2019 0.8390 [0.8198; 0.8587] 13.8 9.6
## Harput 2018 0.8735 [0.8537; 0.8937] 14.2 9.6
## Engelen-VanMelick 2017 1.0333 [1.0011; 1.0666] 7.4 9.5
## Alt 2022 0.9013 [0.8823; 0.9208] 16.3 9.6
## Alt 2022 0.9008 [0.8820; 0.9200] 16.8 9.6
## Gauthier 2022 0.8700 [0.8496; 0.8909] 13.2 9.6
## Mesnard 2022 0.7636 [0.7326; 0.7958] 4.3 9.3
## Raoul 2018 0.6000 [0.4991; 0.7212] 0.2 5.3
## Heijne 2015 0.9150 [0.8753; 0.9565] 3.8 9.3
## Heijne 2015 0.9550 [0.9157; 0.9960] 4.2 9.3
## measure_3
## Joseph 2020 Slow Isokinetic <= 1 year
## Raoul 2019 Slow Isokinetic <= 1 year
## Harput 2018 Slow Isokinetic <= 1 year
## Engelen-VanMelick 2017 Slow Isokinetic > 1 year
## Alt 2022 Slow Isokinetic > 1 year
## Alt 2022 Fast Isokinetic
## Gauthier 2022 Slow Isokinetic <= 1 year
## Mesnard 2022 Slow Isokinetic <= 1 year
## Raoul 2018 Slow Isokinetic <= 1 year
## Heijne 2015 Slow Isokinetic > 1 year
## Heijne 2015 Fast Isokinetic
##
## Number of studies combined: k = 11
##
## ROM 95%-CI z p-value
## Common effect model 0.8878 [0.8802; 0.8955] -27.08 < 0.0001
## Random effects model 0.8727 [0.8197; 0.9291] -4.26 < 0.0001
##
## Quantifying heterogeneity:
## tau^2 = 0.0105 [0.0054; 0.0555]; tau = 0.1026 [0.0732; 0.2355]
## I^2 = 95.0% [92.8%; 96.6%]; H = 4.49 [3.72; 5.41]
##
## Test of heterogeneity:
## Q d.f. p-value
## 201.44 10 < 0.0001
##
## Results for subgroups (common effect model):
## k ROM 95%-CI Q I^2
## measure_3 = Slow Isokinetic <= 1 year 6 0.8534 [0.8433; 0.8637] 54.35 90.8%
## measure_3 = Slow Isokinetic > 1 year 3 0.9371 [0.9218; 0.9527] 50.50 96.0%
## measure_3 = Fast Isokinetic 2 0.9114 [0.8944; 0.9287] 5.94 83.2%
##
## Test for subgroup differences (common effect model):
## Q d.f. p-value
## Between groups 90.66 2 < 0.0001
## Within groups 110.78 8 < 0.0001
##
## Results for subgroups (random effects model):
## k ROM 95%-CI tau^2 tau
## measure_3 = Slow Isokinetic <= 1 year 6 0.8168 [0.7493; 0.8904] 0.0105 0.1026
## measure_3 = Slow Isokinetic > 1 year 3 0.9481 [0.8705; 1.0327] 0.0054 0.0735
## measure_3 = Fast Isokinetic 2 0.9248 [0.8736; 0.9790] 0.0014 0.0377
##
## Test for subgroup differences (random effects model):
## Q d.f. p-value
## Between groups 7.12 2 0.0284
##
## Details on meta-analytical method:
## - Inverse variance method
## - Restricted maximum-likelihood estimator for tau^2
## - Q-Profile method for confidence interval of tau^2 and tauforest_cc_function(hs_ecc_meta, overall = FALSE, prediction.subgroup = TRUE, weight.study = "same", hetstat = FALSE)
Between Person
Only 3 studies report between-person comparisons for eccentric, across quite varied timepoints. No meta-analysis run. Plot here to show effects.
8. Bivariate redo
# Split based on contraction types
joined <- quad_joined %>%
filter(!(cohort == "Li 2022" & measure == "quad isk con 30"),
!(cohort == "Siney 2010" & measure == "quad isk con 90")) %>%
filter(!(cohort == "Bailey 2019" & measure == "quad isk con 300"),
!(cohort == "Drocco 2017" & measure == "quad isk con 240"),
!(cohort == "Kyritsis 2016" & measure == "quad isk con 300"),
!(cohort == "Laudner 2015" & measure == "quad isk con 300"),
!(cohort == "Pamukoff 2018" & measure == "quad isk con 120"),
!(cohort == "Tourville 2014" & measure == "quad isk con 300"),
!(cohort == "Welling 2020" & measure == "quad isk con 300"),
!(cohort == "Welling 2018b" & measure == "quad isk con 300"),
!(cohort == "Zult 2017" & measure == "quad isk con 120")) %>%
filter(!(cohort == "Karanikas 2005" & measure == "quad mvic 0"),
!(cohort == "Labanca 2018" & measure == "quad mvic 30"),
!(cohort == "Labanca 2016" & measure == "quad mvic 30")) %>%
group_by(study, measure, timepoint) %>%
filter(length(type) == 2) %>% # make sure to only include studies who measure both effects
ungroup() %>%
group_by(study, measure, type) %>%
slice(which.min(abs(timepoint_mean - as.numeric(12)))) %>% # take one effect per study, preferably closest to 12 months if multiple
ungroup() %>%
mutate(measure_2 = factor(measure_2, levels = c("Slow Isokinetic", "Isometric", "Fast Isokinetic"))) %>%
group_by(study, type) %>%
arrange(study, measure) %>%
mutate(row = row_number()) %>%
filter(row == 1) %>%
select(-row) %>%
ungroup() %>%
filter(!is.na(vi)) %>% # remove any missing data
group_by(study) %>%
filter(n() > 1) %>%
ungroup()
joinedV <- vcalc(vi,
cluster = study,
type = type,
rho = 0.7,
data = joined,
checkpd = TRUE)
joined_mv <- rma.mv(yi,
joinedV,
mods = ~ type - 1, # remove intercept to generate estimate for each.
random = list(~ type | study),
struct = c("UN"),
data = joined,
cvvc = "varcov") # need this to use later with matreg
# Split based on contraction types
hs_joined <- hs_joined %>%
filter(!(cohort == "Lee 2019" & measure == "hs isk con 60 - deep"),
!(cohort == "Kim 2011" & measure == "hs isk con 90 hyperflex"),
!(cohort == "Siney 2010" & measure == "hs isk con 90"),
!(cohort == "McRae 2013" & measure == "hs isk con 60 supine"),
!(cohort == "Li 2022" & measure == "hs isk con 30")) %>%
filter(!(cohort == "Drocco 2017" & measure == "hs isk con 240"),
!(cohort == "Kyritsis 2016" & measure == "hs isk con 300"),
!(cohort == "Laudner 2015" & measure == "hs isk con 300"),
!(cohort == "Tourville 2014" & measure == "hs isk con 300"),
!(cohort == "Welling 2020" & measure == "hs isk con 300"),
!(cohort == "Welling 2018b" & measure == "hs isk con 300")) %>%
filter(!(cohort == "Hu 2020" & measure == "hs mvic 30 hhd"),
!(cohort == "Hu 2020" & measure == "hs mvic 60 hhd"),
!(cohort == "Ardern 2010" & measure == "hs mvic 30"),
!(cohort == "Ardern 2010" & measure == "hs mvic 105")) %>%
group_by(study, measure, timepoint) %>%
filter(length(type) == 2) %>% # make sure to only include studies who measure both effects
ungroup() %>%
group_by(study, measure, type) %>%
slice(which.min(abs(timepoint_mean - as.numeric(12)))) %>% # take one effect per study, preferably closest to 12 months if multiple
ungroup() %>%
mutate(measure_2 = factor(measure_2, levels = c("Slow Isokinetic", "Isometric", "Fast Isokinetic"))) %>%
group_by(study, type) %>%
arrange(study, measure) %>%
mutate(row = row_number()) %>%
filter(row == 1) %>%
select(-row) %>%
ungroup() %>%
filter(!is.na(vi)) %>% # remove any missing data
group_by(study) %>%
filter(n() > 1) %>%
ungroup()
hs_joinedV <- vcalc(vi,
cluster = study,
type = type,
rho = 0.7,
data = hs_joined,
checkpd = TRUE)
hs_joined_mv <- rma.mv(yi,
hs_joinedV,
mods = ~ type - 1, # remove intercept to generate estimate for each.
random = list(~ type | study),
struct = c("UN"),
data = hs_joined,
cvvc = "varcov") # need this to use later with matregQuadriceps - Bivariate
summary(joined_mv)##
## Multivariate Meta-Analysis Model (k = 54; method: REML)
##
## logLik Deviance AIC BIC AICc
## 34.2578 -68.5157 -58.5157 -48.7594 -57.2113
##
## Variance Components:
##
## outer factor: study (nlvls = 27)
## inner factor: type (nlvls = 2)
##
## estim sqrt k.lvl fixed level
## tau^2.1 0.0276 0.1662 27 no casecontrol
## tau^2.2 0.0188 0.1372 27 no within
##
## rho.cscn rho.wthn cscn wthn
## casecontrol 1 - 27
## within 0.7920 1 no -
##
## Test for Residual Heterogeneity:
## QE(df = 52) = 1452.7544, p-val < .0001
##
## Test of Moderators (coefficients 1:2):
## QM(df = 2) = 48.1098, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## typecasecontrol -0.2131 0.0335 -6.3614 <.0001 -0.2787 -0.1474 ***
## typewithin -0.1790 0.0268 -6.6891 <.0001 -0.2315 -0.1266 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1corplot_function(joined_mv)## [[1]]
##
## estimate se zval pval ci.lb ci.ub
## intrcpt 0.3077 0.0898 3.4260 0.0006 0.1317 0.4837 ***
## casecontrol 0.6538 0.1111 5.8827 <.0001 0.4360 0.8717 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## [[2]]
Hamstring - Bivariate
summary(hs_joined_mv)##
## Multivariate Meta-Analysis Model (k = 38; method: REML)
##
## logLik Deviance AIC BIC AICc
## 35.8751 -71.7502 -61.7502 -53.8326 -59.7502
##
## Variance Components:
##
## outer factor: study (nlvls = 19)
## inner factor: type (nlvls = 2)
##
## estim sqrt k.lvl fixed level
## tau^2.1 0.0279 0.1671 19 no casecontrol
## tau^2.2 0.0030 0.0552 19 no within
##
## rho.cscn rho.wthn cscn wthn
## casecontrol 1 - 19
## within 0.6631 1 no -
##
## Test for Residual Heterogeneity:
## QE(df = 36) = 449.4270, p-val < .0001
##
## Test of Moderators (coefficients 1:2):
## QM(df = 2) = 28.0679, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## typecasecontrol -0.0643 0.0400 -1.6060 0.1083 -0.1427 0.0142
## typewithin -0.0652 0.0135 -4.8395 <.0001 -0.0915 -0.0388 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1corplot_function(hs_joined_mv)## [[1]]
##
## estimate se zval pval ci.lb ci.ub
## intrcpt 0.7316 0.0649 11.2662 <.0001 0.6043 0.8588 ***
## casecontrol 0.2190 0.0692 3.1625 0.0016 0.0833 0.3547 **
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## [[2]]